Statistical Theory of Asymmetric Damage Segregation in Clonal Cell Populations

Asymmetric damage segregation (ADS) is ubiquitous among unicellular organisms: After a mother cell divides, its two daughter cells receive sometimes slightly, sometimes strongly different fractions of damaged proteins accumulated in the mother cell. Previous studies demonstrated that ADS provides a selective advantage over symmetrically dividing cells by rejuvenating and perpetuating the population as a whole. In this work we focus on the statistical properties of damage in individual lineages and the overall damage distributions in growing populations for a variety of ADS models with different rules governing damage accumulation, segregation, and the lifetime dependence on damage. We show that for a large class of deterministic ADS rules the trajectories of damage along the lineages are chaotic, and the distributions of damage in cells born at a given time asymptotically becomes fractal. By exploiting the analogy of linear ADS models with the Iterated Function Systems known in chaos theory, we derive the Frobenius-Perron equation for the stationary damage density distribution and analytically compute the damage distribution moments and fractal dimensions. We also investigate nonlinear and stochastic variants of ADS models and show the robustness of the salient features of the damage distributions.


Introduction
Many unicellular organisms among bacteria and yeasts proliferate by binary fission so that each mother cells divides into two seemingly identical daughter cells. However, more careful examination of cell lineages [1,2] revealed that the progeny is in fact slightly asymmetric, one daughter cell has a slightly longer lifespan than the other. This difference, at least in rode-shaped bacteria E. coli, is apparently rooted in the fact that every bacterium is itself slightly asymmetric since its two poles always have a different age: the old pole that existed in its mother cell and the new pole that was created when its mother cell divided. The "age" of a cell can be defined as the age of the old pole (the number of generations it has been in existence). It was found that the daughter cell that inherits the older pole from its mother, grows slower and divides slightly later than the daughter cell that inherits the newer pole. A plausible hypothesis that can explain such a correlation is that the cells gradually accumulate damaged proteins which aggregates near their poles, and therefore the daughter cell that inherits the older pole, will also inherit a larger fraction of damaged proteins hitherto accumulated in the mother cell. Direct visualization of protein aggregates in growing cell lineages corroborates this conjecture [3][4][5][6].
A number of computational and theoretical studies addressed the dynamics of asymmetric damage segregation (ADS) in a growing microbial population and their implications for the overall population fitness [7][8][9][10][11][12][13][14][15]. In [9,11] it was demonstrated numerically that asymmetry accelerates the average growth of the population as a whole. In [14], a kinetic description of damage segregation was developed on the basis of a transport equation for the time-dependent damage distribution function that was applied both to deterministic and stochastic damage synthesis and segregation. Using that equation, moments and the average population growth rate were computed analytically in the limit of small asymmetry. However, a comprehensive understanding of the structure of damage distribution in a population of asymmetrically dividing cells is still lacking. In this paper, we focus on deterministic asymmetric division and analyze this structure using a general Frobenius-Perron equation for the at-birth damage distribution function. It turns out that if the rules controlling damage accumulation and inheritance are deterministic, the system governing the damage distribution is analogous, and for certain class of linear damage accumulation and division rules exactly equivalent, to the Iterated Function System (IFS) known in the theory of fractals [16]. Exploiting this analogy, we show that the asymptotic stationary distribution of damage is indeed fractal and find the spectrum of its fractal dimensions. For more general nonlinear models of damage accumulations and segregation we analyze the structure of the damage distribution functions numerically and show that the its fractal nature is robust. We discuss also deterministic irregularity of the damages in a lineage and relate it, for the linear damage redistribution rule, to chaotic properties of the IFS. If the damage accumulation and/or segregation have stochastic components, the distribution smears out but remains multi-peaked.

Deterministic Models of ADS
In this section we introduce a class of models for the asymmetric damage segregation that will form the basis of our theory. We will suppose that a cell is created at time t 0 and divides at time t 0 + τ . We denote the instantaneous cell damage (it can be, for example, the amount of damaged proteins, which is assumed to be a real number) D(t), where t 0 ≤ t ≤ t 0 + τ . The initial cell damage (inherited from the mother cell) is x = D(t 0 ), and the final damage just before division is y = D(t 0 + τ ), where τ is the cell lifetime (the latter is constant in some setups or damage-dependent in other formulations).
There are three components in each model that determines the distribution of damage within the population of growing and dividing cells: 1. Damage gain specifies temporal dynamics of damage D(t; x) in a cell starting from its value x at birth at t = t 0 to division at t = t 0 + τ . We will denote the final cell damage as y = D(t 0 + τ ; x).

Examples:
a. Vedel et. al model [11]. In this model it it assumed that every cell adds a fixed amount of damage λ over its lifetime, y = x + λ. The damage is added with the constant rate γ, so that for a lifetime τ this rate is γ = λ/τ . The evolution of the damage is thus D(t; x) = x + λ(t − t 0 )/τ .
b. Damage obeys a linear differential equationḊ = βD + γ, where the term βD describes auto-catalytic production or degradation of damage (depending on the sign of β), with the initial condition D(t 0 ) = x.
2. Lifetime τ might be constant (damage-independent), or it can be a functional of the damage D(t). Since in deterministic models at any moment of time the current damage D(t; x) is pre-determined by the initial value x at the instant of birth t 0 , we can assume that the lifetime is a function of the initial damage, τ (x).

Examples:
a. In model [11], the lifetime is assumed to be a linear function of the cell damage at the division y. Because y is determined by the damage at birth x, one can generally write τ (x) = τ 0 + µy(x). b. Chao model [9]. Here it is assumed that the lifetime is the time when some product P , whose synthesis rate depends on the damage,Ṗ = 1−sD(t; x) with P (t 0 ) = 0, reaches a threshold value P 0 . This leads to a nontrivial dependence of the cell lifetime on its initial damage τ (x), see Eq. (37) below.
3. Damage inheritance: A deterministic rule according to which the damage of a mother cell is distributed between two daughters (no damage losses),

Examples:
a. In model [9], a linear mapping x 1,2 = 1±a 2 y with a constant asymmetry parameter 0 < a < 1 is adopted. b. In model [11], a nonlinear mapping x 1,2 = 1±a(y) 2 y with the asymmetry parameter being a function of the damage of the mother cell, a(y), is suggested.
All model variants listed above are deterministic, however we will generalize our description to account for stochasticity in ADS in Sec. 8.

Qualitative description of the deterministic cell population dynamics
It is straightforward to simulate any specific ADS model numerically, starting, for example, from a single cell with a certain initial damage. Fig. 1a,b illustrate these dynamics for a model with the fixed damage gain (rule 1a of the previous Section), linear lifetime function (2a), and linear damage inheritance rule (3a). In the following we will label this model as 1a2a3a. Fig. 1a corresponds to the constant lifetime (µ = 0) and Fig. 1b to the damage-dependent lifetime (µ > 0). Note that since the number of cells in such a simulation grows exponentially with time, it cannot be extended very far. As an alternative, one can simulate a Moran-type process with a constant population of N cells when at each cell division, another random cell is removed from the population, so the population size remains constant. The statistical features of both processes are the same in the "thermodynamic limit" N → ∞, where one can neglect finite-size fluctuations. We can also expect them to be similar for finite N at the time when the size of a growing population is also N .
The main difference between the cases in panels (a) and (b) is that for the damage-independent lifetime all cells divide simultaneously, while for damagedependent lifetimes, synchrony of division times is quickly lost. In both cases, after a short transient, the distribution of damage in just-born cells appears to reach a broad and multi-peaked stationary state (Fig. 1c,d) with a finite support between certain x min and x max . It is straightforward to find x min and x max for a given model as the asymptotic limit of inheriting only smaller or only larger fractions of mother's cell damage within a lineage. For example, for Model 1a2a3a, x min = λ 1−a 1+a , x max = λ 1+a 1−a . It is, however, not trivial to interpret these simulations. Indeed, the process is fully deterministic: an initial damage of a cell fully determines damages and fission times of all its descendants. But it is not a classical deterministic dynamical system because the number of cells grows, and we cannot represent the time evolution of the damage values in all descendant cells as a single trajectory. However, while each mother cell has two daughters, the daughter cell damage uniquely specifies the damage of its mother, at least for the model 1a2a3a, if parameters a and λ are also known. The reason for this is that the damages of daughters receiving larger and smaller fractions of the mother's damage belong to non-overlapping sets (x min , λ) and (λ, x max ) (this is not the case for all ADS models, see below). Thus, for the model 1a2a3a we can construct a deterministic 'back-in-time' map for the damage x(n) of the just-born cells in n-th generation and (e) with (f) we observe that the process quickly to the asymptotic long-time regime. A detailed characterization of the early transients is beyond the scope of this paper. x(n-1) x(n) ( Fig. 2a) and generate a unique "lineage" damage trajectory back in time, as illustrated by Fig. 2b where we concatenated the damage time courses in different generations. The main conclusion of Fig. 2b is that the deterministic nature of the model notwithstanding, the time course of the damage is highly irregular, which suggests that this process may be approached from the viewpoint of the chaos theory. The mapping (1) is indeed an expanding piecewise linear map which is a prototypical model of deterministic chaos [17]. This immediately explains the irregularity of the lineage damage trajectory. In fact, relation (1) can be considered a Poincaré map for the continuous trajectories of the full time-dependent damage dynamics. The latter is described not by differential equations, but rather by a combination of a continuous evolution of damage during the cell lifetime interrupted by a discrete transition from mother to daughter by the rule of damage inheritance. In this sense our hybrid continuous-discrete system resembles an integrate-and-fire model for firing neurons, where a continuous voltage increase (current integration) is combined with a jump rule (firing of a spike and a reset of the voltage).
We stress here that although a back-in-time trajectory can be always extracted from a forward-in-time simulation of a population of dividing cells, not in every case this back-in-time trajectory can be associated with a deterministic dynamical system like the map (1). It requires the uniqueness of the back-intime map, i.e. for any initial damage of a daughter cell there should be only one possible value of the mother's initial damage. A detailed analysis (Sec. 5) shows that this is not always the case. For example, for model 1b2a3a with β > 0, two branches of the back-in-time map similar to Fig. 2b) overlap, and so neither forward-in-time, nor back-in-time lineage trajectories can be computed by iterating a deterministic 1D map.
The essential ingredient of any ADS model is the dependence of cell lifetime on its damage. Continuing with our analogy between the mother-to-daughter damage inheritance and a Poincaré map, one can interpret the lifetimes as Poincaré return times. The properties of these times affect the regularity of the damage trajectories Fig. 2a; the corresponding notion in the chaos theory is phase coherence [18]. One may distinguish a generic case of damage-dependent lifetimes from a degenerate case of fixed lifetimes. For fixed lifetimes, the time evolution of the damage is only partially irregular: Divisions occur at regular time intervals (Fig. 1a), and there is no mixing for the continuous-time dynamical process. All lineages starting from the same cell undergo divisions simultaneously, and so the mean population damage continues to oscillate indefinitely (see Fig. 3a, red line). In other words, all the cells remain "in phase". In contrast, for a generic damage dependence of lifetimes, the intervals between sequential divisions are chaotic (since the initial damages of cells in a lineage are chaotic), and the continuous-time process is mixing (Fig. 1b). Thus, different lineages starting from the came cell decorrelate (in other words, the phases of different cells become scattered) and the mean population damage eventually settles into a steady state (see Fig. 3a, blue line). Another way to quantify this difference is to compute the normalized autocorrelation function of the damage trajectory where the angular brackets denote averaging over different lineages 1 .
In Fig. 3b we compare this autocorrelation function for model 1a2a3a with µ = 0 and µ = 0. One can see that while for damage-dependent lifetimes (µ = 0) this function decays to zero, for constant lifetimes it initially decays (reflecting the irregularity of the damages in just-born cells), but asymptotically at large time lags it oscillates periodically without decay, reflecting phase coherence of the corresponding damage trajectory. For the latter case, the autocorrelation function can be computed analytically (see Appendix A).

Kinetic description of the cell population dynamics
In the previous section we demonstrated that the dynamics of the asymmetric damage segregation in cell lineages is irregular, thus it appears appropriate to describe them statistically. In this section, we formulate a kinetic description of ADS in the thermodynamic limit of a large number of cells. It is similar to the recently developed transport approach [14,15] but differs in some details.
As we have argued above, in many aspects the damage dynamics is analogous to the dynamics of chaotic oscillators, and we will follow this analogy in the construction of the statistical theory. It is convenient to introduce a phase of the cell cycle φ that changes linearly from φ = 0 (birth) to φ = 1 (division) with the rate ω (which is simply the inverse of the cell lifetime τ ) that generally may depend on the initial cell damage x: Next, instead of D(t; x) above, we can introduce the damage during the cell lifetime as a deterministic function of the initial state x and the phase φ, D = F (x, φ). This function is related to the function D(t; x) above as Let us now introduce a two-dimensional population density N (x, φ, t), so that N (x, φ, t)dxdφ is the number of cells in a small interval of the initial damages dx and in a small interval of phases dφ. With a proper normalization, the quantity N (x, φ, t) can be interpreted as a probability density in the space of deterministic variables x, φ. We stress here that the irregularity of the dynamics is not needed for this: even regular deterministic processes (e.g. populations of periodic oscillators) can be described via a probability density in the phase space, with an additional assumption of random initial phases in an ensemble of copies [21]. This quantity is similar to the distribution ψ(u, x, t) introduced in [14], which instead of our phase φ depends on the actual cell age u measured in units of time. We believe that our distribution of cell ages in terms of the phase φ (i.e. the fraction of the cell's lifetime) has its advantages, especially when the lifetime varies strongly from cell to cell.
Instead of dealing with an exponentially growing population of dividing cells, it is convenient to introduce cell death, so that a certain fraction of cells ΓN (x, φ, t)dt dies within a small time interval dt irrespective of x or φ. The death rate Γ in principle can be arbitrary, but if we choose it exactly balance the (yet unknown) average asymptotic proliferation rate, the total number of cells will remain constant, and the distribution N (x, φ, t) should eventually approach a stationary state N (x, φ). It is easy to see that the death rate Γ found from the stationarity condition of the asymptotic solution for a model with death is exactly the growth rate of a proliferating cell population without death. Alternatively, we can assume a Moran-type process [22] when at the time of every binary fission, one other randomly chosen cell is taken away, so the total number of cells remains precisely constant at all times, not just asymptotically and in the thermodynamic limit. It can be shown (see Appendix B) that asymptotically in the large-time limit the kinetic equation for the Moran-type process is equivalent to the one presented here.
The transport equation for N (x, φ, t) in the system with cell death follows from the probability conservation: where the second term on the l.h.s. describes cells advection from phase 0 to phase 1 with phase velocity ω(x) according to (3), and the r.h.s. term describes cell death. This equation formally coincides with the transport equation for the normalized damage distribution in an exponentially growing population [14] with Γ playing the role of the population growth rate instead of the death rate in our model. We assume that the death rate Γ is indeed chosen properly to balance the population growth. In this case, and if the initial phases in a population are random, the probability distribution asymptotically approaches a stationary state described by the equation 2 The solution of this equation describes the exponentially decaying flux of cells from φ = 0 towards φ = 1: To close the system and to ensure stationarity, we need to add the boundary (transition) condition balancing the incoming flux of the daughter cells N (x, 0)ω(x) at φ = 0 with the outgoing flux of mother cells N (x, 1)ω(x) at φ = 1 (due to the cell division, the number of cells leaving the interval at φ = 1 has to be exactly the half of the number of cells entering at φ = 0, thus factor 2 in Eq. (14) below) . The connection is determined by the rule specifying how the damage at the end of the life cycle y = F (x, 1) is redistributed within daughter cells. For a general damage inheritance rule (see 3. in section 2) we havē For monotonous functions F and f , depending on their specific forms, a daughter cell with initial damagex can descend from one of two mothers with initial damages . The balance condition for the fluxes thus can be written as Dividing both sides by dx we arrive at the basic Frobenius-Perron operator for the damage distribution in a population of dividing cells at the beginning of their cell cycle Alternatively, we can write this Frobenius-Perron equation for the product where g 1,2 (x) = f 1,2 (F (x, 1)) and we dropped the bar over x. It is easy to see that properly normalized P (x) is the density of just-born cells (the fraction of the whole population of cells that that are born within a small unit time interval with initial damage x). This equation can also be deduced directly from the general equation for the stationary initial damage distribution function derived in [14] P where f (x|x ) is the probability to have initial damage x in a cell if the initial damage in its mother cell is x . Eq. (12) can describe both deterministic and stochastic damage distribution scenarios. To obtain our deterministic (11) one just needs to take g(x|x

The conservation of the number of cells requires
This self-consistency condition together with Eq. (11) uniquely determines the death/growth rate Γ.
Asymptotically, the solution of the Frobenius-Perron equation (10) (provided the death rate Γ satisfies (13)) converges to the stationary damage distribution at the moments of birth N (x, 0). It is easy to see that this distribution has a finite support between x min and x max which are determined from equations Let us introduce the total number of cells (here we use also (7)): Then the stationary probability distribution density over initial damages and phases is We can use this two-dimensional density to obtain the density for the distribution of the damages (now not initial one, but just observed at all possible phases), by the following expression Here we used expression (4) which relates the damage at the phase φ to the initial damage x. We can also obtain he distribution over the phases w(φ) as the marginal distribution by integrating the two-dimensional density (16)

Linear damage accumulation, constant lifetime
In this section we use the kinetic theory developed in the previous section to analyze the statistical properties of the model with a fixed lifetime τ = T and linear mappings with constant A and B where for definiteness, we assume A < B.
It is easy to see that linear transformations (18) with A + B = 1 correspond to rules 1a and 3a of Sec. 2, where The same linear mappings also can be deduced from rule 1b, however in this case, integration of the ODE for the damage within the cell cycle yields and so for non-zero β, A + B = 1. If β < 0, then A + B < 1, and if β > 0, then A + B > 1. Note that for the damage to remain bounded at all times it is required that The fixed lifetime may correspond either to rule 2a with µ = 0 (then T = τ 0 ) or to rule 2b with s = 0 (then T = P 0 ).

Fractal properties of the damage distribution
The discrete dynamics of x according to two linear transformations (18) is a well-known mathematical object called Iterated Function System (IFS), see [23][24][25]. The particular case A = B is often called Bernoulli convolution [26]. IFS is the simplest mathematical model generating fractals [16]. To illustrate this, we first write down the equation for P (x) that follows from the general Frobenius-Perron equation (11): where we also took into account the condition e ΓT = 2 which immediately follows from (14) for Numerical iteration of the operator (19) for arbitrary A and B generically yields a fractal distribution, as exemplified by Fig. 4a where we used the same parameters as in Fig.1a 3 . This distribution is virtually indistinguishable from the distribution obtained in direct simulations of the underlying ADS model (Fig.1c).
The fractal properties of this distribution can be summarized as follows (see [25] for the current state of the theory).
Case without overlap A + B ≤ 1.. This is the simplest case where there is no overlap of two linear branches in the mapping (18). In this case one can characterize the invariant measure by generalized fractal dimensions d q (see, e.g., [17]; traditionally for these dimensions capital letters are used, but in this paper this notation is reserved for the damage). A standard scaling argument [17] leads to the exact parametric expression Most important are the box-counting dimension d 0 and the information dimension d 1 . One can see that when the damage of just-born cells is conserved A + B = 1, the box-counting dimension is one, i.e. the support of the measure is the full interval [x min , x max ]. In other words, there are no voids and the set of possible damages is not a classical Cantor set. In contrast, if the initial  damage is partially dissolved, i.e. A + B < 1, then d 0 < 1 and the set of possible damages is a Cantor set. In both cases the information dimension d 1 < 1 (with the exception of a trivial degenerate situation of symmetric segregation A = B = 1/2, when the measure is uniform).
Case with overlap A + B > 1.. In this case the two branches of (18) overlap. This situation has long been a conundrum for the measure theory, and only recently it has been partially clarified [23][24][25]. In particular, Ref. [25] analyzed generalized dimensions in the range 0 ≤ q ≤ 1. For typical AB > 1/4, all these dimensions are d q = 1, what means that the distribution is continuous with a finite density P (x). For AB < 1/4 and A + B > 1, there is a "phase transition" in dependence on q: there is a critical value q * , beyond which the expression (20) holds, while below this value D q is a fractional-linear function of q, with d 0 = 1.
A more detailed understanding exists for the symmetric case A = B, which is called Bernoulli convolution. Here, it has been proven that for almost all values of A > 1/2, the invariant measure is absolutely continuous [27]. At exceptional points (so-called Pisot numbers) the distribution is fractal, but the information dimension is very close to one (see recent estimates in [28,29]). We refer to Ref. [30] for a nice illustration of the densities for different values of A > 1/2.

Singularities at the ends of the interval
One can see from Fig. 4 that the behavior of the damage distribution density is very different at the opposite ends of the interval, i.e. close to x min and x max . While there is a sharp peak at x min , density at x max nearly vanishes. To understand this structure, let us consider a vicinity of a fixed point, for definiteness the fixed point x min . Let us denotex = x − x min . Then the first branch is g 1 (x) = Ax, while the second branch Thus in the vicinity ofx = 0, the branch g 2 takes the iterate of a smallx far away from the vicinity of the fixed point, while the branch g 1 brings the iterate of a smallx even closer to zero according to the stable linear mappingx → Ax. Correspondingly, the Frobenius-Perron equation (19) for the invariant density around the left fixed point reduces to since the second term on the r.h.s. of (19) gives zero contribution nearx = 0. Seeking the solution in the form P (x) = Cx γ−1 we obtain The power γ defines the singularity of the density P (x). One can see that the critical value of the mapping slope is A = 1/2. If A < 1/2, then γ < 1 and the distribution has a diverging peak near x min . If A > 1/2, then γ > 1, and the density vanishes at the fixed point. The same condition holds for the right fixed point at x = x max : here the density has a peak for B < 1/2, otherwise the density vanishes. For the case depicted in Fig. 4, we have A < 1/2 and B > 1/2, and thus the density has a peak near x min and vanishes near x max . The second mapping branch g 2 (x) "transports" the boundary peak at x min to the position g 2 (x min ), which is subsequently split into two more which are located at g 1 (g 2 (x min )) and g 2 (g 2 (x min )) and so on, so that an infinite set of peaks appears (with progressively smaller amplitudes since at each new iteration every peak is split into two). It is interesting to note that the condition for the absence of peaks on both ends of the interval is that both A, B > 1/2, which coincides with the above-mentioned necessary condition for a continuous density AB > 1/4.

Moments and autocorrelations of the damage distribution
It is remarkable that although the distribution of damages in just-born cells is fractal, its statistical characterization in terms of moments is quite simple and can be computed analytically for arbitrarily large asymmetry a. This, of course, is due to the linearity of the basic equation (19).
We can compute moments M k of the distribution of P (x) or the initial damage N (x, 0) by multiplying Eq.(19) by x k and integrating. After straightforward algebra, we get for the first two moments Note that these averages over a distribution of damage in a population at a given time are different from the lineage distributions mentioned in Sec. 3 although the two are related [19,20].
Using the linearity of the governing kinetic equation (19) we can also compute the normalized auto-correlation function of the damages in just-born cells of different generations in a single lineage, where the argument is the integer generation number. Calculation of the C(m+ 1) with relation between x(n + m + 1) and x(n + m) given by (18), leads to a recursion This yields exponentially decaying correlations In fact, the full continuous-time correlation function presented Fig 3 can also be calculated analytically (see Appendix A).

Phase-averaged distribution of damage
Because of the linear relation P (x) = T −1 N (x, 0), we can use the properly normalized P (x) (i.e. P (x) dx = 1) in the two-dimensional damage density that in this case is factorized: Thus, the marginal distribution over the phases is exponential We can also compute the phase-averaged distribution W (D) from Eq.(17), however the result depends on the explicit form of the function D = F (x, φ).
Linear damage growth.. Substitution of F (x, φ) = x + λφ in (17) yields the following expression for the distribution of the damage integrated over the phases As seen in Fig. 4, unlike P (x), W (D) is continuous. Using this expression, one can explicitly calculate the moments of the distribution damages: we get the following expression for the phase-integrated distribution of the damage With β → 0 and λ = T γ this reduces to the expression (25).

Linear damage accumulation, damage-dependent lifetimes
Here, we consider the same linear model with the mother-daughter damage inheritance relations (18), but assume that cell lifetimes depend on the damage, τ (x).

Stationary damage distribution
We again start with the general Frobenius-Perron equation (11) and rewrite it for the linear damage distribution model (18): or Unlike the previous section, the value of Γ here is unknown and needs to be determined from the conservation of the total probability. Numerically, it can be implemented iteratively in two different ways. In the first, we solve Eq. (28) iteratively starting from an arbitrary initial distribution P 0 (x), but at each iteration, we find P k (x, Γ) for a set of values Γ. The normalization condition P k (x, Γ) dx = 1 can be considered as equation for Γ, the root of which is determined numerically. Thus, the proper normalization is ensured at each iteration. As a result of these iterations of the Frobenius-Perron operator (27), we obtain a sequence of densities and of values of Γ, both of which converge. The corresponding limit are the stationary density and the corresponding stationary death rate. The second method is to start with some initial guesses for P (x) and Γ, compute new P (x) using (27), then compute S = P (x, Γ) dx (which is generally not 1) and update Γ = Γ · S. Then compute P (x) in the next iteration, and do it until both P (x) and Γ converge to their asymptotic values. The second method is more computationally efficient and precise, however its convergence generally is not guaranteed.
The stationary distribution in Fig. 4b obtained by solving the FP equation agrees very well with direct numerical simulation shown in Fig. 1d. Qualitatively, it is also fractal, however, since formally equation (27) does not correspond to a classical IFS, we cannot rely on the corresponding mathematical theory and compute fractal dimensions of the invariant measure analytically. Nevertheless, one can evaluate the generalized dimensions numerically. We illustrate this in Fig. 4d, where we show the spectrum of fractal dimensions of asymptotic damage distributions for several values of parameter µ (while other parameters remain fixed).

Cumulant expansion of the damage distribution
In this section we present an approximate analysis of the Frobenius-Perron equation (28) for the case of the linear dependence of lifetimes on damage, τ (x) = T (1 + µ(x − x 0 )). To simplify the calculation, without loss of generality we assume that x 0 = x µ=0 , the average damage of new-born cells for a fixed lifetime (see expression (22)).
The main idea of the analysis below is to explore cases of weak lifetime variability across the population. As the expression for τ (x) shows, this occurs if µ(x − x 0 ) is small in the range of x values of the whole population, i.e. when either the parameter µ is small, or the deviations (x − x 0 ) are small. A good measure of these deviations is the variance (23) that is proportional to (A−B) 2 . Thus, the asymptotic analysis presented in this Section is valid either for weak dependence of lifetimes on the damage (small µ) or weak asymmetry of damage segregation (small ε = (A − B) 2 ).
Below we will only sketch the theory, see Appendix C for a full derivation. The method is based on expanding the characteristic function Using the Frobenius-Perron equation (27), we can easily write the equation for C(k): Taking the logarithm of both sides and introducing the cumulant-generating function F (k) = log C(k), we obtain the following equation (30) Next we substitute the general cumulant expansion and arrive at (32) Equating terms at powers of k, we obtain a system of equations for the cumulants.
Let us first briefly discuss the case of constant lifetimes µ = 0. In this case, equations in the order m contain only cumulants with indices m ≤ m. Thus, the cumulants can be calculated sequentially starting from c 1 . In fact, this procedure is equivalent to the ad hoc derivation of moments in Sec. 5.3. Unfortunately, this property is lost for µ = 0. However, as we shall argue below, in many interesting cases one can perform a truncation of the infinite system of equations for the cumulants. Below we use a three-cumulants truncation: setting all c m , m > 3 in (32) to zero gives a system of four equations for unknown Γ, c 1 , c 2 , c 3 : Inspection of these equations reveals that there are indeed two potentially small parameters, justifying the truncation: 1. Small non-isochronicity of lifetimes, i.e. small parameter µ. Cumulants in this case do not need to be small. As one can conclude from (33), ΓT is in fact represented by a power series in µ where higher orders in µ come with higher cumulants. Also the higher cumulants enter (34) multiplied with powers of µ. This allows for calculating Γ approximately, as a series in µ. This case was also treated in [14] by a direct moment expansion of their general transport equation.

Small higher cumulants. This occurs if the difference |x
For small µ, we calculated the death rate Γ up to O(µ 3 ): ΓT = ln2 + µ 2 Γ 2 + µ 3 Γ 3 (the term ∼ µ is absent due to the proper choice of the central value x 0 = x µ=0 in the definition of τ (x)). The explicit expressions for Γ 2,3 are given in Appendix C. In Fig. 5a,b we compare this approximate analytical expression with numerics. For small ε, the three-cumulants truncation only determine Γ up to O(ε), because c 4 ∼ ε 2 is missing in (33). The calculated approximate expression ΓT = ln2 + Γ 1 ε (see Appendix C for the explicit expression for Γ 1 ) is compared with numerics in Fig. 5c,d. The comparison presented in Fig. 5 illustrates an excellent agreement between numerics and the asymptotic theory. The results obtained here are also in agreement with earlier numerical and analytical findings [9,11,14,31]. Specifically, while the mean damage x is a decreasing function of both µ and ε, the trends for the distribution variance are opposite. The increase in segregation asymmetry ε obviously leads to a wider distribution of damages (greater variance). As we have seen in the previous section, for µ = 0 the mean damage of a population is independent of ε. A non-zero µ gives selective advantage to cells with less damage and thus reduces the mean for finite compared with ε = 0. The increase in µ for a fixed ε gives greater selective advantage to cells with smaller damage and therefore reduces both the mean and the variance of the distribution. Also, the population-averaged growth rate Γ is an increasing function of both µ and ε, that independently contribute to the spread among cells' lifetimes, at least when these two parameters are small, i.e. for a small spread of lifetimes in a population. Note that for more complex models relating cell growth and damage [14], the dependence of the population growth rate on the asymmetry may be non-monotonous and switch from beneficial to detrimental at some finite ε.

Nonlinear deterministic ADS models
Above we focused on models where many relevant statistical properties can be found analytically, either explicitly or as perturbative expressions for some small parameters. The main simplifying assumption is a linear redistribution of damages, which corresponds to linear Iterated Function Systems. For more generic nonlinear models we do not have an analytical theory, and the goal of this section is to show that numerical simulations reveal features similar to those in the analytically tractable cases.

The Chao model
As a representative example here we consider the Chao model [9]: • Damage gain: The instantaneous damage in every cell grows linearly with time, D(t; x) = x + γ(t − t 0 ) • Lifetime is the time when some product P whose synthesis is suppressed by the damage according toṖ = 1 − sD(t; x) with P (t 0 ) = 0, reaches a certain threshold value P 0 . One can easily see that the lifetime τ is a solution of the quadratic equation (1 − sx)τ − sγ 2 τ 2 = P 0 . Note, that this model of lifetime assumes that the product sD remains small so thaṫ P > 0.
Noteworthy, although the damage inheritance itself is linear, the effective map of damages of just-born cells is nonlinear because of the nontrivial lifetime dependence on damage. Below, for compatibility with the theory above, we use A = (1 − a)/2, B = (1 + a)/2. The dependence τ (x) comes as a solution of the quadratic equation above Thus, the damage just prior to division is related to the initial damage as a nonlinear transformation and thus the two branches of the mapping of damage from one generation to the next are g 1 (x) = Ay(x), g 2 (x) = By(x).
Minimal and maximal values of x are determined from x min = Ay(x min ) and x max = By(x max ): Figure 6 illustrates the transformation g 1,2 (x) for A = 0.35, B = 0.65, s = 0.4, P 0 = 0.3, γ = 0.85. As before, it consists of two branches, but now the branches are nonlinear, furthermore, for the given parameters they overlap also along the vertical axis (this means that the mother's initial damage is not uniquely determined by the daughter's initial damage). Simulations of this model yield the distribution of x shown in Fig. 7a. As already discussed in Section 5.1, maps like Fig. 6 with a large overlap result in a peaky but formally non-fractal (at least in the sense of absence of voids) distribution. However, the existing literature about fractal properties of IFS is mainly restricted to the linear case, thus exact statements about the fractal properties of the distribution are hardly possible. One can see from Fig. 7a that the minima of the distribution are separated from zero, which is an indication of a density without voids. In this figure we also present the autocorrelation function to confirm irregularity of the damage time dependence.   Here we discuss some general statistical properties of distributions produced by nonlinear IFSs, without connecting them directly to particular damage production and segregation models. IFS are characterized by two functions g 1,2 (x), so that in every iteration a given x produces two new two states x 1,2 = g 1,2 (x), thus the number of states x in each generation k grows exponentially as 2 k . We assume that g 1,2 (x) do not intersect, and g 1 (x) < g 2 (x). The interval of possible values of x is limited by the two fixed point where g 1,2 (x) intersect the diagonal: x min = g 1 (x min ) and x max = g 2 (x max ). Without loss of generality we assume that x is normalized such that x min = 0 and x max = 1, i.e. the distribution of x lies within a unit interval 0 ≤ x ≤ 1.
To illustrate the variety of distributions that can be generated by nonlinear IFSs, we consider parabolic functions g 1,2 (x): : [here parameters A, B define the slopes of f 1 (x) and f 2 (x) at the fixed points x = 0 and x = 1 and parameters A 2 , B 2 characterize nonlinearity of the functions].
In Fig. 8 we present six typical regimes, for which we numerically computed the distribution densities (P (x) and the corresponding cumulative distributions defined as W (x) = x 0 P (x )dx . We also compute the autocorrelation functions of the iterative sequences of x. In all cases the autocorrelation function 4 decays exponentially, so that the damage level along a fixed lineage has strong chaotic properties.
Fractal properties of the damage distribution strongly depend on the presence of an overlap or a gap between the ranges of values of functions g 1,2 on the interval 0 ≤ x ≤ 1: these ranges are 0 ≤ g 1 (x) ≤ g 1max and g 2min ≤ g 2 (x) ≤ 1. Three situations without a gap are depicted in panels (a,b,c). The case of panel (a), with a vanishing gap, is qualitatively similar to the standard linear IFS with parameters A + B = 1, discussed in Section 5.1 above. The measure is fractal, but it does not have voids. Panels (b,c) of Fig. 8 show two situations with an overlap, they differ by the stability properties of the fixed points at x min , x max . Since near the fixed point nonlinear mappings can be linearized, we can used the results of Section 5.1 (Eq. (21)): if A > 1/2 and B > 1/2 (case of panel (b)), the density at both fixed points vanishes. This results in a rather smooth distribution which strongly resembles a Gaussian hump. On the contrary, for A < 1/2, B < 1/2 (case of panel (c)), there are singularities at the fixed points, which are "transported" along the interval by functions g 1,2 , so that the final structure contains a sequence of peaks. We are not aware of any statements/conjectures about absolute continuity of the measure in such a nonlinear case.
Three situations with a gap g 1max < g 2min are depicted in panels (d-f). Cases (d,f) resemble a classical Cantor set with a large void in the center and an hierarchy of smaller voids in the density distribution (these voids correspond to the horizontal intervals of cumulative distributions in Fig. 8). In panel (f) the fractal is symmetric (following the symmetry of the functions g 1,2 ), while the distribution in panel (d) has a peak at x min and a vanishing density at x max , again according to the values A, B. In contradistinction to these cases, for non-monotonous functions g 1,2 in panel (e) we observe just one large void. The distribution density in this case is concentrated in two separated regions close to zero and close to one, and inside these regions there are no additional voids (there the distribution is presumably continuous).
Returning to the asymmetric damage distribution problem, let us note that the borderline no gap/no overlap IFS of the linear model 1a2a3a mostly analyzed in the earlier sections is a degenerate situation caused by the assumption that the difference between the damages in two sister cells stays the same during the cells lifetimes. In general, even though during the cell division the damage is conserved, f 2 (y) = 1 − f 1 (y), different damage gains in two sister cells (possibly caused by the autocatalytic nature of damage production or the lifetime dependence on the damage) lead to either a gap or an overlap between the two branches. When the difference between the end damages in the sister cells is smaller than their initial difference, (y 2 − y 1 < x 2 − x 1 ), the corresponding IFSs have an overlap and typically produce a continuous, in some cases even rather smooth distributions. If, on the other hand, the difference in damages becomes stronger over the lifetime (y 2 − y 1 > x 2 − x 1 ) there is a gap between the values of the two branches, and the distribution is typically fractal with voids (and even with a hierarchy of voids in the case of a Cantor-type measure).

Stochastic damage distribution and segregation
In biology, all processes are stochastic, due to extrinsic environmental fluctuations and intrinsic randomness of biochemical reactions that is particularly important on a sub-cellular scale. Randomness can in principle affect all three rules that constitute a model of asymmetric damage segregation. The fraction of the mother's final damage y that a daughter cell inherits can be stochastic and governed by the conditional distribution w 1 (x|y). Note that since for a single mother cells there are two daughter cells, it is convenient to normalize this distribution as w 1 (x|y)dx = 2. Furthermore, because for two daughters' damages satisfy x 1 + x 2 = y, the conditional distribution should be symmetric w 1 (y − x|y) = w 1 (x|y).
Damage gain and lifetime can also fluctuate and together they are specified by a two-dimensional distribution w 2 (y, τ |x) conditioned on the initial damage x (note that y and τ are generally not independent, since longer-lived cells may on average accumulate more damage). However, it appears plausible to assume that segregation is statistically independent of the damage accumulation. Thus, the stochastic Frobenius-Perron equation for the probability distribution P (x) can be written as follows where again the normalization condition dx P (x ) = 1 yields the equation for the growth rate Γ. This equation is a more general form of the self-consistent equation for the distribution P (x) derived in [14] [they postulated a deterministic relationship between x and τ ]. Note that if we take the probability distributions in the form w 1 (x|y ) = δ(x−f 1 (y ))+δ(x−f 2 (y )); w 2 (y , τ |x ) = δ(y −F (x ))δ(τ −τ (x )) , we recover the Frobenius-Perron equation (11) for the fully deterministic case.
If the stochasticity only affects the damage inheritance while the damage accumulation and the lifetime are deterministic functions of the initial damage, y(x), τ (x), we can substitute w 2 (y , τ |x ) = δ(y − y(x ))δ(τ − τ (x )) in (39) and arrive at the Fredholm integral equation of the second kind that is equivalent to the self-consistent equation of [14]. If the probability distribution w 1 (x|y ) is continuous but still close to the two-peaked form, for example with small spread σ (w 0 is the normalization constant to satisfy w 1 (x|y )dx = 2), the damage distribution that was fractal in noise-free system, becomes continuous but still highly irregular (see Fig. 9). On the other hand, damage distributions that were continuous and smooth in deterministic limit, are much more robust agains noise (data not shown). Note that the distribution of damage in all (and not in just-born) cells at a given time is smooth even even when the distribution of initial damages is fractal, and so it is also very robust agains noise. The distribution of stochastic damage segregation can also be single-peaked. As an example, let us consider a truncated Gaussian and as before assume linear functions for the damage gain y(x) = x + λ and lifetime τ (x) = τ 0 + µ(x + λ). Figure 10 shows distributions of initial damage  P (x) for two different magnitudes of randomness σ. For small σ, the distribution P (x) is narrow and appears close to a Gaussian. In fact, it is easy to see that for small σ, when the truncation in (42) can be neglected, the solution of the FP equation (40) is Gaussian. Substituting in (40) it is easy to solve for x 0 , ∆, and Γ (see Appendix D). Fig. 11 shows the mean and standard deviation of damage and the population growth rate as functions of σ obtained from this approximate solution and directly from numerical simulations of the underlying stochastic ADS model. In agreement with earlier fundings [14,31], the population-averaged growth rate Γ increases with the randomness, which should not be very surprising since inheritance randomness effectively creates asymmetry in damage separation. In fact, asymmetry in cell division is not necessary for creating and exploiting the selective advantages of a broad damage distribution. Even if divisions are deterministic and symmetric, w 1 (x|y ) = 2δ(x − y /2), but the damage accumulation is stochastic, the damage distribution will have a finite width and the average growth rate will also be greater than in a purely deterministic symmetric case. In this case, the Frobenius-Perron equation (39) simplifies to Let us consider a simple illustrative example where the damage accumulation is described by a truncated Gaussian distribution and the lifetime τ is a deterministic linear function of the final damage y, τ = τ 0 + µy: For small σ we again can ignore the truncation and solve the Frobenius-Perron equation (44) analytically by substituting the solution in the same Gaussian form (43) (see Appendix D). The resulting x 0 , ∆, and Γ as functions of σ are shown in Fig. 12 superimposed with the results of direct numerical simulations of the underlying stochastic process.

Discussion
Aging of microbial populations has been a subject of active research in recent years [32]. While superficially fissioning bacterial cells appear immortal and divide symmetrically via binary fission, more close inspection reveals a slight phenotypic asymmetry that has been attributed to the asymmetric inheritance of damaged and aggregated proteins accumulated in the mother cell, among its daughters [4]. Importantly, a daughter cell inheriting a greater fraction of ancestor's damage, replicates slower than its sibling that inherits a smaller fraction of the damage. Cells that have a long line of ancestors inheriting lesser fraction of the damage, become "rejuvenated" and divide more often thus producing more offspring than "age" cells that have predominantly ancestors with greater fractions of damage accumulate larger amounts of damage. Asymmetry in damage inheritance has also been found in yeast [3,6] and even in higher eukaryotes [33]. A number of conceptually simple models have been proposed to describe this branching process [7][8][9][10][11][12][13][14][15]. Simulations and analysis of these models revealed that asymmetric damage segregation might be evolutionally preferred because it accelerates the mean population growth by letting rejuvenated cells, however most of these results were found from direct simulations.
We re-examined these models focusing on understanding the mathematical and statistical properties of the damage distribution in populations of asymmetrically dividing microbes. In a very broad class of models that encompasses both deterministic and stochastic ADS rules, the asymptotic damage distribution in the beginning of cell cycle can be described by a Frobenius-Perron-type equation (39) where the rules of damage accumulation and inheritance are encoded in the transition probabilities w 1 (x|y), w 2 (y, τ |x) for damage inheritance and damage accumulation, respectively. For deterministic transition rules, the damage distributions are broad and highly irregular resembling a fractal set. The mappings of the initial damage from generation to generation are equivalent to to the Iterated Function Systems, and for linear mappings the fractal dimensions of the stationary damage distributions can be computed analytically. Stochasticity in damage accumulation and segregation smoothes out the fine fractal structure of the distributions, however for small noise, the distributions remain highly irregular, and their moments as well as the average growth rates of the population remain robust. We expect that experiments with fluorescent labeling of damaged proteins will reveal a complex multi-peaked structure of their distributions predicted by our theory. The overall width of the distribution will give us a quantitative measure of the asymmetry in damage inheritance, while the presence and magnitude of the distribution peaks (fractality) will characterize the degree of deterministic vs. stochastic asymmetry in damage distribution. Another characteristics which is potentially measurable in experiments and has attracted less attention in previous literature, is the autocorrelation function of lineages (2). Its structure sheds light on mixing properties of the ADS and can be compared with theoretical expressions, Here we would like briefly discuss the applicability of our stationary solutions for finite-size and growing populations. For a theoretical treatment, it is convenient to balance exponential population growth with cell death so the population reaches a stationary state in a statistical sense thus allowing for the customary characterization of stationary statistical processes (invariant distribution density, correlation function, Lyapunov exponent, etc.). It is easy to see that population-normalized damage distributions in exponentially growing and stationary populations are described by identical Frobenius-Perron (or transport) equations with death rate Γ playing the role of mean growth rate Λ (see [14]). In finite-size population starting from a single cell (like in Fig. 1 (a, b)) or a small group of cells, additional dependencies on the initial state and on the system size appear. We did not address them in the present study, but these transient aspects (which appear to be relevant for experimental observations) definitely deserve future study.
Phenotypic variability in clonal populations is often invoked as a useful bethedging strategy that improves the population survival chances in adverse environmental conditions [34]. Phenotypic variability comes in many forms and can be caused by many factors such as mutations, multistability in underlying gene networks, noise in gene expression, or post-transcriptional processes, etc. While not any phenotypic variability is beneficial, the asymmetric damage segregation that maintains low damage in a fraction of the population to offset continuous accumulation of damaged proteins in all cells has been shown to help colony survival in stressful conditions when cells with sufficiently high of damage become mortal [35,36].
(D.3) Equating O(x 2 ) terms in the exponents on both sides, we obtain ∆ 2 = (∆ 2 + 4σ 2 )/4 or ∆ 2 = 4σ 2 /3 . These dependences are shown by blue solid lines in Fig. 11. Similarly, for the case of symmetric partition and Gaussian damage accumulation distribution, we substitute (45) in (44) and again assume σ is small and ignore the truncation: These dependences are shown by blue solid lines in Fig. 12. Interestingly, formulas (D.10), (D.11) for the mean initial dmage x 0 and the growth rate Γ coincide with expressions (42),(D.6) while the formula (D.9) for the width of the initial damage distribution ∆ is different from (D.4).